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ABSTRACT 

A comparison of numerical methods utilized by the finite 
element technique for solving a nonlinear nuclear reactor 
dynamics problem was conducted. Using the Crank-Nicolson, 
DVOGER (Gear) and Implicit Gear methods, the results showed 
the Implicit to be the superior method investigated. This 
is based on the fact that all three methods yielded the same 
steady state solutions; but, the Implicit Gear method used 
significantly less CPU time and comparable storage to Crank- 
Nicolson. This was particularly apparent as the degrees of 
freedom were increased. In addition, the transient solution 
in all cases was better than that obtained in Crank-Nicolson 
and compared favorably to that of Gear’s method. 

The other noteworthy result was in the effect of the 

error criterion on solution. It was shown that for a range 
- 4 

of error from 10 to 1.0, the steady state solution value 
remained the same. This results in a significant reduction 
in computer processing time since the time required decreases 
substantially as the error conditions imposed are relaxed. 
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I. INTRODUCTION 



A. PURPOSE 

This research project has been undertaken to compare 
several numerical methods of solving a nonlinear nuclear 
reactor dynamics problem. Three methods have been investi- 
gated in this thesis, including the Crank-Nicolson, the 
DVOGER (Gear) and the Implicit Gear methods of solution. 

A nuclear reactor dynamics problem with temperature- 
dependent feedback, when it entails either a non-homogeneous 
or multi-region reactor, results in a nonlinear field equa- 
tion in space and time. This problem does not lend itself 
to solution by analytical means [5, 6]. However, when the 
physical and neutronic properties of the problem are known, 
a model can be formulated using the finite element method 
which will yield the transient and steady state flux solu- 
tions. In particular, the partial differential equations 
investigated were of the form 

a -||- = bV^ ijj + cip - o)ip ^ (1) 

where a, b, c, and co are constants and ij; (r, z, t) is the 
flux. The finite element method reduces Equation (1) to 
the system of ordinary differential equations 
N . N 

j=l A ij *j (t) = j=l B i j *j (t) 1 = N W 

when the nonlinear term is linearized. Nguyen and Salinas [5] 
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and more recently Olsen [6] discuss the problem and methods 
of solution. 

In this work, a comparison of three numerical integration 
schemes for Equation (2) was made. The comparisons will be 
made on several bases: computer storage requirements, com- 

puter processing time, rapidity with which solution was 
obtained and the relative accuracy of the solution. 

At steady state, it is expected that all three methods 
will provide the same solution value. This is due to the 
fact that at steady state the time derivative of the flux 
( iJO , Equation (2), is zero. 

In order to explore the relative value of the various 
methods, the model has been discretized into finite element 
grids of various degrees of freedom (DOF) . The effect on 
the solution by finer discretizations of the finite element 
has been investigated to determine if the various methods of 
solution are similarly influenced to provide a better solution 
for a finer mesh. To have a method of solution relatively 
independent of mesh size would greatly reduce computer storage 
and processing time requirements since a larger grid with 
fewer elements could be utilized. 

In order to test the flexibility of the various equation 
solvers, an initial disturbance was introduced at different 
points in the model. This provided both a check of the 
ability of the method to accept a random disturbance as well 
as information on how rapidly a solution is obtained with a 
particular disturbance input. Two nodal points of the system 
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were considered, a point at the origin and a point on the 
core-reflector interface. This was done to determine if the 
tracking ability was consistent throughout the model. 

Additional comparisons investigated include the effects 
of the convergence criterion on the solution for all methods 
and the effects of the size of the time step on the Crank- 
Nicolson method. 

Modifications have been made to the programs provided in 
Ref. [6] which was the basis of this research project. In 
that work, core properties were arbitrarily assigned to the 
reflector elements at the interface. This occurred because 
the material and nuclear properties were provided at the nodal 
points rather than by elements. The properties were intro- 
duced in this manner as a means of computer storage reduction. 
Regardless of the mesh size, there were always fewer nodes 
than elements. In this project, all properties were provided 
on an element basis to eliminate this discrepancy while at 
the same time sacrificing the minor storage savings it repre- 
sented. 
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I I . DATA GENERATORS 



A. GENERAL 

The data generators formulated in this work are useable 
only for regular rectangular grids (Figures 1-3) . Having 
selected the number of horizontal and vertical points for 
the discretization of the model, the data generators will 
provide the numbering of each nodal point, the horizontal 
and vertical position of each node and the nodal neighbors 
of each node. 

In addition, the outer boundary nodes, at which the 
dynamic flux is zero, will be numbered last in order to 
reduce the number of equations required by the particular 
method of solution being used. This will substantially 
decrease the storage requirements. For example, in a one 
hundred thirty-two node discretization, only one hundred ten 
equations will be solved. 

B. PROPERTY INPUTS 

1 . Purpose 

A simple data generator has been provided which will 
produce a data deck containing the physical properties of 
the reactor core and reflector in the format required by the 
Crank -Ni co Is on and DVOGER (Gear) methods. These properties 
are neutron velocity, neutron diffusion length, neutron 
absorption cross-section, neutron fission cross-section and 
reactivity temperature coefficient. 



14 



This generator will accommodate any size mesh and 
will process an indefinite number of grids simultaneously. 

The user must provide this routine with a data deck 
which contains the total number of elements in the grid and 
the type of each element, core or reflector. This has been 
simplified by the designation of all core elements by ITYPE=0 
and reflector elements by ITYPE^l. 

2 . Programming Notes 

a. The property values in the two regions must be 
provided as an integral part of the program. 

b. The type of element, core or reflector, must be 
provided. 

c. The routine will process an unlimited number of 
models. It, therefore, must be halted when all data has been 
utilized. This is done by specifying the final value of the 
number of elements in an IF-STOP statement. 

3 . Parameters 

a. NEL - number of elements in the discretized model 

b. ITYPE - type of element; ITYPE=0 is a core element 
and ITYPE=1 is a reflector element. 

c. D - vector containing the diffusion length of the 
core and reflector elements. 

d. SGA - vector containing the absorption cross- 
section of the core and reflector elements. 

e. SGF - vector containing the fission cross-section 
of the core and reflector elements. 

f. V - vector containing the neutron velocity. 
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g. ALPHA - vector containing the reactivity- 
temperature coefficient. 

C. NODAL POINT COORDINATES AND ELEMENT 
NODAL POINT CONNECTIVITY 

1 . Purpose 

Nodal point positioning in the model is readily 
obtained in data deck form from this routine. Once the model 
to be investigated has been discretized into the finite ele- 
ment grid desired, the user provides the number of vertical 
and horizontal nodal points and their dimensional position 
along the axes. By the use of a nested loop, the nodal points 
are consecutively numbered and assigned the proper coordinate 
dimensions . 

The element connectivity, that is, the nodal points 
for each element, is also resolved by the use of a nested 
loop and several counters. Each iteration will yield two 
elements with their respective nodal point boundaries. This 
will continue until the nodal points for each element have 
been computed. 

The output of this routine will consist of two data 
decks. The first will provide the nodal point with its ver- 
tical and horizontal coordinates. The second will yield the 
element nodal point connectivity. 

Omitted from the element connectivity data deck, but 
required by the thesis program, is the type of element, core 
or reflector. This must be added to the deck individually. 
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Example of output (See Figure 1) : 



nodal point 
1 

element 

1 



nodal point coordinates 
horizontal position vertical position 
0.0 0.0 



element connectivity 
nodal point connectivity 
19 10 



type element* 
0 



*Note: must be added manually to data card. 



2 . Programming Notes 

a. The maximum number of vertical and horizontal 
nodal points in the discretized model must be provided. 

b. The horizontal and vertical positions of the 
discretization must be provided. 

c. The boundary points of the model will be numbered 
such that these nodal points are the last in the numerical 
sequence . 

d. The routine is written such that it will process 
an indefinite number of models. Therefore, it must be halted 
when all the data has been utilized. This will be accomplished 
by specifying the maximum number of vertical points in the 
model in an IF-STOP statement. 

3 . Parameters 

a. NH - maximum number of horizontal nodal points 
in the discretized model. 

b. NV - maximum number of vertical nodal points in 
the discretized model. 
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c. X - vector containing the horizontal positions 
of the NH points. 

d. Y - vector containing the vertical positions of 
the NV points. 

e. SYSNOD - vector providing the total number of 
nodal points. 

f. R - vector providing the radial position of the 
nodal points. 

g. Z - vector providing the vertical position of 
the nodal points. 

h. ELNOD - array providing the nodal point boundaries 
for each element. 

i. NEL - the total number of elements in the discre- 
tized model. 

D. NODAL NEIGHBOR CONNECTIVITY 
1 . Purpose 

This routine produces a data deck for the Crank- 
Nicolson and Implicit Gear methods containing each nodal 
point and the nodal points connected to it by the finite 
element discretization. 

The regular rectangular grids are such that no more 
than six nodal points contribute to any one; therefore, an 
Nx7 array can be formulated for nodal neighbor connectivity. 
For example, nodal point 26 in the 112 element grid would be 
stored as follows (Figure 4) : 

26 17 18 27 35 34 25 



18 



and nodal point 58 



58 49 50 59 57 0 0. 

As written, zeros are inputted whenever a nodal point does 
not have six neighbors to maintain the symmetry of the matrix. 
This could, however, be compacted further by reading the data 
as a single vector with a starting point counter to indicate 
when the next nodal point has been reached in the sequence. 
This would eliminate the need for the zeros to be supplied. 

This routine will provide a data deck for any 
rectangular grid and for an indefinite number of grids. 

2 . Algorithm 

This algorithm provides a dense, compacted matrix 
which reduces the storage size required by the main program. 
The use of a finite element method allows a certain amount of 
compacting, i.e., the banded matrix. This is a function of 
the finite element size, that is, the greater the number of 
nodal points, the larger the band. The size of the band is 
determined by the largest difference between nodal neighbors 
plus one. For example, in Figure 1, node 11 has a maximum 
difference of 20 - 2 = 18, while at node 16 the difference 
is 43 - 7 = 36 which is the largest difference in the forty- 
five node system. As the number of nodal points increases, 
this maximum band width also increases. In Figure 2, the 
maximum difference is at node 16 (67 - 7 = 60] ; and, in 
Figure 3, the maximum difference of 124 - 10 = 114 exists 
at node 22. Therefore, although the size of the system is 
reduced from NxN to Nxq, the size is not constant and may 
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not be small. In general, if one used banded storage, the 
numbering would proceed consecutively in the direction of 
fewest nodes. However, this would require the identification 
of those nodes on the boundary since they are of constant 
value and do not require integration. 

In the particular scheme of discretization utilized 
in this project, no nodal point has more than six nodal point 
contributors. This allows the matrix to be compacted further 
from Nxq to Nx7. 

3 . Programming Notes 

a. The total number of nodal points, the maximum 
number of nodal point contributors and the maximum number of 
horizontal and vertical nodal points must be provided. 

b. The routine will satisfy any rectangular grid. 
However, it must be instructed when a sufficient number of 
nodal points have been generated. This is accomplished by 
using the total number of nodal points desired in an IF-GO TO 
statement . 

c. This routine requires a positive means of stopping. 
This is accomplished by the use of the maximum number of 
vertical nodal points in the last data set in an IF-STOP 
statement. 

d. The output of this routine will place the central 
nodal point first with the contributors following. 

4 . Parameters 

a. NV - number of vertical nodal points in the model. 

b. NH - number of horizontal nodal points in the model. 
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c. NVH 



total number of nodal points in the model. 



NVH = NV x NH, 

d. LCON - maximum number of contributing nodal points. 
(LCON = 7) . 

e. MNOD - an array, NVH x LCON, providing the central 
and contributing nodal points. 
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III. CRANK-NICOLSON METHOD OF SOLUTION 



A. DESCRIPTION 

The Crank-Nicolson formulation is a numerical method 
originally presented by J. Crank and P. Nicolson in 1947. 
Crank-Nicolson, being an implicit method, does not require 
the inverse of the matrix to be calculated. Therefore, 
advantage is taken of the sparse matrix inherently provided 
by the finite element method. If the inverse of the sparse 
matrix is formed, it will be full. 

The Crank-Nicolson method will solve the set of linear 
differential equations represented by 

jii A ij *j (t) = j?i B ij *j (t > i=1 > (3) 

Crank-Nicolson, after some algebra, yields 



+ f - At 


- r 


+ - At 


At 


Lj 


2 



The implicit system (4) may be solved by an iterative process; 
in this case, the Gauss-Seidel method is used to solve the 
set of simultaneous algebraic equations (4) formulated. Of 
all the stable numerical methods in the case of single step 
implicit, Crank-Nicolson has the smallest truncation error [2] . 
The Gauss-Seidel method of iteration is continuously using the 
newest solution values allowing convergence to the solution 
to be the most rapid. 
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B. EFFECT OF DEGREES OF FREEDOM 



As might be expected, as the grid size became finer, the 
computer core requirements increase. In addition, more time 
is required for the problem to reach a steady state solution. 
At the same time, as shown in Figures 5-7, the solution values 
converged, as the mesh size became finer, to a more accurate 
solution. 

As shown in Table I and Figure 8, storage and processing 
time, that is, CPU time to reach a steady state solution, in- 
creased markedly between a forty- five node grid and a one 
hundred thirty-two node grid. This yielded a fifty-nine 
percent increase in storage and better than a five hundred 
percent increase in processing time. At the same time there 
is an improvement in the solution of only four and one-half 
percent . 

Comparison to a seventy-two node grid yielded far more 
satisfying results. At seventy-two nodes, there is an eighty 
percent increase in CPU time and a sixteen percent increase 
in storage requirements while obtaining a two and one-half 
percent improvement in the solution. 

Similar results were obtained in problem time to solution. 
For one hundred thirty- two nodes, the time more than doubles; 
while for seventy- two nodes, the time was increased by twenty 
percent . 

C. EFFECT OF ERROR CRITERION 

The error criterion used throughout this thesis informs 
the integration routine when it has achieved sufficient 
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agreement between iterates (i.e., when it can stop and go 
to the next iteration step) . 

In Crank-Nicolson, the error criterion simply consisted 
of the difference in successive Gauss-Seidel iterates divided 
by the current iterative value. 

The error criterion was varied from 10" 1 to 10~^ to 
observe the effects on the solution and computer requirements. 
For Crank-Nicolson the steady state solution value remained 
the same regardless of the error criterion. This demonstrated 
that the extra time required to satisfy a tighter error cri- 
terion at each iteration was not necessary to reach a satis- 
factory steady state solution. As was expected, the computer 
processing time did increase significantly as the error 

criterion was made smaller. An unusual result occurred when 

- 4 

the error criterion was set to 10 . This particular value 

resulted in a processing time less than that required for 
10 \ It would appear that this might have been the result 
of two contributing factors. With the closeness of the error, 
each time an iteration was performed, it was using a better 
solution, and each new time step also had a better solution 
value. Although there were many more time steps required 
for this error criterion, the steady state value was rapidly 
approached because fewer iterations were required at each new 

time. These results are shown in Table II. These results 

- 1 - 4 

for an error criterion between 10 and 10 do not imply 
that this would be the result obtained for all problems of 
this type. 
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D. EFFECT OF INCREASE IN TIME STEP SIZE 

In this method, the time step can be increased or not 
depending on the solution. If a solution is not obtained 
with a particular time step, that same time step is reduced; 
and the routine attempts to attain solution with the smaller 
time step. This continues until either a satisfactory solu- 
tion value is obtained or the built in default value is 
reached which terminates the routine. When a solution is 
attained, the routine compares the number of iterations 
required to reach that value to that of the previous time 
step and to a specified number of iterations. If the current 
number of iterations is less than either of those numbers, 
the program increases the time step by an input constant 
value . 

The change made to the initial time step as the steady 
state was approached was the critical value if the solution 
was to converge to a steady state. It appeared that Crank - 
Nicolson was particularly sensitive to the amount the time 
step was increased as the solution was approached (Figure 9) . 
The method was unable to recover if, when nearing the knee of 
the curve, the size of the increase was such that the steady 
state value was exceeded. Once the solution was passed, the 
results indicated a diverging oscillation about the steady 
state value. 

Several values of the initial time step were utilized in 
all grids to attempt to locate the steady state. The results 
of these trials are provided in Table III. It is noteworthy 
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that if the increase was greater than 1.2 times the initial 
step size, the steady state would never be achieved. Whether 
this would be the case for all grids is not obvious. This is 
particularly evident in view of the results obtained by 
Olsen [6], who utilized an increase of 1.5 times the initial 
step for a thirty-eight node grid. It would, therefore, seem 
that a trial and error approach would be necessary until an 
increase in step size resulted in a steady state solution. 



26 



IV. DVOGER (GEAR) METHOD OF SOLUTION 



A. DESCRIPTION 

The DVOGER (Gear) routine is an IMSL library routine 
which integrates a system of explicit first order differ- 
ential equations. Within this library routine is the capa- 
bility of solving both stiff and nonstiff systems by selection 
of an indicator. In solving the nonstiff system, the Adams 
predictor-corrector is utilized. For the stiff system, Gear’s 
predictor-corrector is used. Gear’s method computes the 
Jacobian for the system of ordinary differential equations 
in order to optimize the time step for each integration. 

B. EFFECT OF DEGREES OF FREEDOM 

Gear's method placed a much greater demand on the computer 
in terms of processing time and storage than the other methods 
considered. This was primarily due to three causes: 1) the 

calculation of the Jacobian, 2 ) the transformation of Equation 
(3) to explicit form, and 3) the initial absolute value of 
the solution. The first two items require a substantial 
increase in both CPU time and computer storage. Several 
variances were implemented in this program to determine if 
better use could be made of the method to make it competitive 
with Crank-Nicolson or the Implicit Gear methods. 

When the system was considered nonstiff, the storage 
requirements decreased since the Jacobian was not calculated; 
but, the processing time increased by one order of magnitude 
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when compared to the stiff system treatment due to the small 
time steps taken to solution. The routine calls for an 
initial absolute solution value of one to be supplied. This, 
however, adversely affects the progress of the problem since 
this value of one is initially compared to values of 10*^. 

This limits the time step taken by the routine initially and 
continues to affect the progress until the steady state solu- 
tion is neared. When an initial value of 10 14 was utilized, 
the processing time was decreased by forty percent. This was 
due to the fact that larger time steps were allowed from the 
outset, since the previous solution value (in this case the 
initial value) was comparable to the value obtained in the 
first calculation. 

At best, when utilizing Gear’s method, the processing 
times were two orders of magnitude greater than that required 
for the Implicit Gear method and twenty times the Crank- 
Nicolson processing times for the 132 DOF system, as shown in 
Table I. Core requirements were also increased significantly, 
particularly at one hundred thirty- two degrees of freedom. 

For this grid, DVOGER (Gear) was approximately three times 
Crank-Nicolson and the Implicit Gear core requirements. 

Gear's method provided the same steady state solution 
values as the other two methods, and the transient curves 
were very similar to those obtained in the Implicit Gear 
method (Figures 10-15) . 



28 



C. EFFECT OF ERROR CRITERION 

In Gear’s method, the time step size is adjusted so that 
the single step error estimate divided by the previous maximum 
solution value is less than the error criterion in the 
Euclidean norm. The single step error estimate is a multiple 
of the difference between the predicted and corrected values 
of the variable. 

The error criterion. Table II, was varied for this method 
as it was for the others. The results were the same; as the 
criterion was tightened, CPU time increased. In this case, 
however, the time became excessive very rapidly, more than an 
hour for 10 ^ and almost four hours for 10 ^ for the model 
with 45 DOF. Similarly, the solution values were the same 
regardless of the criterion utilized. 
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V. IMPLICIT GEAR METHOD OF SOLUTION 



A. DESCRIPTION 

This numerical method is particularly useful for the 
solution of large, sparse systems of implicit, stiff differ- 
ential equations [4]. In contrast to Gear's method, the 
Implicit Gear method eliminates the need of an exact NxN 
Jacobian, but rather requires only an exact Nx7 Jacobian. 

This is due to the fact that Equations (3) are handled directly, 
thereby retaining the sparseness of the matrix. 

Gear's predictor- corrector method and the compact matrix 
makes use of storage because the implicit Equations (3) are 
handled directly. This results in very efficient use of the 
computer both in terms of storage and processing time require- 
ments . 

The user is required to provide a subroutine that evalu- 
ates the system of equations being investigated , as well as, 
for efficiency, a subroutine to evaluate the Jacobian of the 
system of equations. 

This program is a modification by Franke [4] to the 
routine DFASUB [7]. One of the major changes to the routine 
is in the treatment of the error. In the Implicit Gear 
method, instead of using the Euclidean norm, the root mean 
square norm is utilized. This is no more than the Euclidean 
norm divided by the square root of the number of components. 

In addition, the maximum value of the component is updated 
before the norm of the relative error is computed. 
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B. EFFECT OF DEGREES OF FREEDOM 

The steady state values obtained in this method were the 
same as those obtained in Crank-Nicolson and DVOGER (Gear) 
as shown in Figures 11-18. This result does not imply that 
these methods will always give identical steady state solu- 
tions. Certainly they do give different transient solutions. 
The major advantages of this method are 1) the computer 
processing time was reduced as the grid size became finer 
and 2) the storage requirements are slightly greater than 
that required in the Crank-Nicolson method due to the compu- 
tation of the Jacobian and the storage of up to seven past 
solutions which control the size of the time step. Implicit 
Gear requires about Nx25 more storage locations than Crank- 
Nicolson. This amounts to about only thirteen thousand bytes 
for the one hundred thirty-two degrees of freedom system 
(single precision) . The magnitudes of the increase in pro- 
cessing time dropped dramatically in this method. For the 
same initial disturbance, the time increased eighty percent 
for the one hundred thirty- two node grid and twenty- four 
percent for the seventy-two degrees of freedom. In addition, 
Implicit Gear gave a more accurate transient solution. The 
results of this method are shown in Table I and comparison 
to the other methods will be made in Chapter VI. 

C. EFFECT OF ERROR CRITERION 

In Implicit Gear, the error criterion is defined as in 
DVOGER (Gear) except that the root mean square of the Euclidean 
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norm and the updated maximum solution value are used in the 
computation. 

The error criterion was varied substantially to determine 

what effect it had on computer storage and time requirements, 

as well as its effect on the accuracy of solution. As shown 

in Table II and Figures 19-20, a wide range of convergence, 

-4 

1.0 to 10 , was investigated with interesting results. 

Expectedly, the processing time did increase as a tighter 
error criterion was required. However, the criterion had 
little effect on the track through the transient solution for 
values of less than 10 When using an error of one-half 

and one, the transient solution varied substantially. The 
same steady state solution was obtained, although at a later 
point in problem time as shown in Figure 19. For the stiff er 

error criterion requirements, the processing time necessary 

- 4 

is doubled by utilizing a criterion of 10 instead of a value 
greater than 10" 1 . This, certainly, did not appear to warrant 
the closer tolerance level. 
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VI . RESULTS 



A. COMPARISON OF TIME AND STORAGE REQUIREMENTS 

With all three methods yielding the same steady state 
solution (Figures 10-18), the comparison of methods became 
one of time and storage requirements placed on the computer 
rather than one based on which provided the best solution. 

In the transient stage. Gear tracked slightly better than 
Implicit Gear, but both Gear and Implicit Gear tracked better 
than Crank-Nicolson (See Figures 10-18). For the three systems 
of equations utilized, the Implicit Gear method performed 
significantly better than either Crank-Nicolson or DVOGER 
(Gear) in CPU time and was comparable to Crank-Nicolson and 
superior to Gear in storage requirements (Table I) . This 
becomes more significant as the number of degrees of freedom 
increase. The Implicit Gear is less sensitive to the increase 
in size of the system of equations as shown in Figure 8. 

Comparing the forty-five and the one hundred thirty-two 
degrees of freedom, the core requirements increased slightly 
and the time doubled for the Implicit Gear method. In con- 
trast, for Crank-Nicolson, the time increased by six times 
and the core by 100K bytes; and, for DVOGER (Gear), there 
was a ten times and 500K byte increase in time and storage 
respectively . 
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B. DEGREE OF FREEDOM EFFECT ON SOLUTION 

Varying the degrees of freedom resulted in a better 
solution value as shown in Figures 5-7 and 21-26. There was 
a four and one-half percent better solution obtained for the 
one hundred thirty-two degrees of freedom compared to the 
forty-five degrees of freedom. This better solution is very 
costly both in terms of computer core requirements and CPU 
time (Figure 8) . 

C. ERROR CRITERION 

Error criterion was varied for all methods to determine 
the effect on solution and computer. As might be expected, 
the computer processing time increased for all methods as 
more rigid tolerance was imposed (Table II) . However, although 
the criterion was made very close, there was no appreciable 
effect on the transient solution until the error criterion 
was greater than 10 1 and no effect on the steady state 
solution. This is significant in that, as stated above, the 
processing time increases greatly as the criterion is tight- 
ened. This shows that for the particular problem considered 
here, some close error criteria yield no advantage, only the 
disadvantage of requiring more time in the computer. 

D. INITIAL DISTURBANCES 

The input of various initial disturbances (central, skewed 
at the core- ref lector interface and uniform throughout the 
core) yielded the same steady state solution value. The 
transient solution varied substantially due to the input 
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position of the disturbance. In most cases, due to the 
magnitude and scope of the initial disturbance, the uniform 
disturbance required the least processing time to reach the 
steady state value. 

E. COMPARISON OF SOLUTION TRACKING 

Two nodal points were investigated graphically for each 
disturbance, integration method and DOF, as shown in Figures 
5-7, 10-18 and 21-26. At these test points, the methods 
performed similarly in all cases in both the transient and 
steady state solutions. 
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VII. TEST PROBLEMS 



A. PHYSICAL MODEL 

A cylindrical reactor with the dimensions and properties 
given in Table IV was used as a model. A radial slice of 
this model was discretized by various finite element grids 
as shown in Figures 1-3. 

B. COMPUTER PROCESSING CONSIDERATIONS 

The programs presented in this thesis. Appendix A, were 
written in the FORTRAN IV language and all computer runs. 

Table V, processed on the IBM 360/67 computer using the 
FORTRAN ' G f compiler. It is expected that the FORTRAN 'H f 
compiler would have supplied results similar to those obtained 
with the FORTRAN ’G' compiler. This was not done because the 
FORTRAN 'H' compiler requires 350K bytes as a minimum core 
requirement. The majority of the runs performed in this 
thesis required significantly less than 300K storage. Single 
precision (six to seven significant digits) was used through- 
out this research. As has been shown [6] , double precision 
solutions are at variance by less than 0.01 percent with 
single precision solutions. 

C. PROBLEM ANALYSIS 

The techniques used by Olsen [6] , modified as described, 
and the Implicit Gear method [4] were utilized to solve the 
problem delineated in Ref. [5] . The problem was approached 
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using finite element discretizations with forty-five, seventy- 
two and one hundred thirty- two degrees of freedom and providing 
an initial disturbance. Three initial disturbances were pro- 
vided as follows: central at the core center (R = 0 cm., 

Z = 0 cm.); uniform throughout the core; and, a skewed dis- 
turbance at the core-reflector interface (R = 60 cm. , Z = 0 
cm. ) . 



D. PROGRAM USAGE 

In order to properly utilize the routines presented in 
this thesis, several points must be considered. 

In the Crank-Nicolson routine, there was difficulty in 
obtaining a steady state solution. The ultimate cause was 
the size with which the previous time step increases. This 
requires a trial and error approach for the particular grid 
size. For this project, an increase of ten or twenty percent 
yields satisfactory results; but, anything larger results in 
a divergent oscillation about the steady state value. 

In all three methods, the dimensioning should be set by 
the particular problem being solved, i.e., determined by the 
degrees of freedom. If this is done, most efficient use will 
be made of the computer, and the user will experience better 
turn around time on the equipment. 

In the use of the Implicit Gear method, the data cards 
for nodal neighbor connectivity must have the contributing 
nodes listed consecutively, with the zeros, used for nodes 
not having six neighbors, being last on the cards. For example, 
in the one hundred thirty-two degrees of freedom system 



37 



(Figure 3) , the nodal neighbor connectivity for node 122 
would be written as 

122 123 11 0 0 00. 

The data deck arrangement for all three methods is given in 
Table VI . 
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VIII. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

It has been shown in this research project that of the 
three methods investigated in the solution of the test problem 
that the Implicit Gear performed in a superior manner in CPU 
time requirements and was comparable to the storage require- 
ments of Crank-Nicolson. In addition, Implicit Gear was the 
least sensitive to the change in the degrees of freedom. 

The method used had no effect on the steady state solution 
value obtained as all three yielded the same results fox the 
same DOF. In the transient solution, Implicit Gear conformed 
very closely to DVOGER (Gear). As shown in Figures 10-15, as 
the DOF increased, the transient solutions moved in the direc- 
tion of the transient solution of Gear's method. 

As the DOF in the discretized finite element was increased, 
a better solution was obtained. This, however, was counter- 
acted by the fact that for the small increase in accuracy 
obtained by the finer mesh, the time and storage need of the 
computer increased significantly. 

A very important result of this project was the effect of 
varying the error requirements of the problem. This yielded 
essentially no change in the accuracy of the transient solu- 
tion for a range of lO** 1 to 10" 4 and no variance in the steady 
state solution for a range from 1.0 to 10 4 . In using these 
methods of integration, the error criterion selected would 
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be a function of the solution desired. If the user is only 
concerned with the steady state solution, the tolerance could 
be relaxed to 1.0. However, if the transient solution is 
needed, an error of 10 ^ appears to be the least rigid value 
which will still provide a satisfactory track to steady state. 

B. RECOMMENDATIONS 

It would be of value to investigate further the Implicit 
Gear method through the use of other problems solving a non- 
linear system of ordinary differential equations. In addition, 
larger degrees of freedom should be attempted to evaluate any 
limitations that may exist in this method of integration. 
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TABLE I 



COMPARISON OF SOLUTION METHODS 



METHOD 


CORE 


(K) 


PROCESSING TIME* 
Central Skewed 


(min. ) 
Uniform 






45 


Nodes 






CRANKO 


158 




1.67 


1.73 


3.57 


DVOGER (GEAR) 


148 




17.8 


19.5 


12.5 


IMPLICIT GEAR 


124 




1.05 


1.00 


1.20 






72 


Nodes 






CRANKO 


184 




3.07 


5.49 


2.38 


DVOGER (GEAR) 


244 




89.5 


114.1 


59.5 


IMPLICIT GEAR 


124 




1.3 


1.30 


1.25 






132 


Nodes 






CRANKO 


252 




9.1 


9.1 


9.1 


DVOGER (GEAR) 


610 




>400 


>400 


>400 


IMPLICIT GEAR 


124 




1.9 


2.0 


2.0 



* Processing time is that required tb reach a steady state 
solution. 
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TABLE II 



EFFECTS 


OF SOLUTION ERROR 


CRITERION 


SOLUTION ERROR 
CRITERION 


45 Nodes 

PROCESSING TIME 
(min. ) 


TIME AT SOLUTION 
(sec.) 


0.1 


CRANK-NICOLSON 

1.8 


3.69 x 10' 5 


0.01 


5.5 


3.70 x 10' 5 


0.001 


14.2 


3.57 x 10" 5 


0.0001 


0.83 


3.87 x 10~ 5 


1.0 


IMPLICIT GEAR 
0.99 


7.56 x 10' 4 


0.5 


0.92 


1.48 x 10' 2 


0.1 


1.05 


7.06 x 10' 5 


0.01 


1.20 


7.64 x 10‘ 5 


0.001 


1.45 


6.36 x 10' 5 


0.0001 


1. 58 


6.44 x 10‘ 5 


0.1 


DVOGER (GEAR) 
18.6 


5.64 x 10' 5 


0.01 


29.. 5 


5.79 x 10‘ 5 


0.001 


78.8 


5.55 x 10' 5 


0.0001 


>300 


- 
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TABLE III 



EFFECTS OF VARYING THE INCREASE OF THE INITIAL TIME STEP 



IN CRANK-NICOLSON METHOD 



INCREASE 


CORE 


PROCESSING TIME 
(min. ) 


SOLUTION TIME 
(sec. ) 


1.1 


158K 


45 NODES 
4.99 


3.18 x 10' 5 


1.2 


1 58K 


1.65 


2.69 x 10' 5 


1.5 


184K 


7.42 


OSCILLATING 


1.1 


2 52K 


72 NODES 
3.5 


4.38 x 10' 5 


1.2 


184K 


3.05 


3.23 x 10' 5 


1.5 


252K 


3.5 


OSCILLATING 


1.4 


252K 


9.1 


OSCILLATING 


1.5 


184K 


6.1 


OSCILLATING 


1.1 


252K 


132 NODES 
9.1 


5.58 x 10~ 5 


1.2 


252K 


8.4 


8.03 x 10' 5 


1.5 


184K 


2.6 


OSCILLATING 
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TABLE IV 



PHYSICAL CONSTANTS 



SYMBOL 



COMPUTER 

SYMBOL 



DEFINITION 



VALUE 



C-N, D / IG 
R R 

R c R 

H c Z 

H Z 

v V/VELOCT 

D D/DSUBF 

c 

D r D/DSUBC 

Z SGA/SGMAAF 

ac 

Z ar SGA/SGMAAC 

v ZNU 

I £c SGF/SGMAF 

Z £r SGF/ - 

e FISFAC 

hA/V HBAR 

a ALPHA 



Total radius 

Core radius 

Core height 

Total height 

Neutron velocity 

Neutron diffusion 
coefficient (core) 

Neutron diffusion 

coefficient 

(reflector) 

Neutron absorption 

cross-section 

(core) 

Neutron absorption 

cross-section 

(reflector) 

Number of neutrons 
per fission 

Neutron fission 

cross-section 

(core) 

Neutron fission 
cross-section 

Fission energy 

Modified convection 
heat transfer 
coefficient 

Reactivity temper- 
ature coefficient 



90 cm. 

60 cm. 

160 cm. 

220 cm . 

4 . 8 x 10 7 cm/sec 
0.913 cm. 

1.20 cm . 

0.01401 cm' 1 
0.008 cm' 1 
2.54 

0.008 cm' 1 
0.0 cm 1 

7.652 x 10' 12 cal/ fis 
0.0632 cal/cm 3 -sec-°C 

1 x 10' 5 / °C 
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RUN 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 



TABLE V 



LIST OF COMPUTER RUNS 



METHOD NODES DISTURBANCE COVERGENCE TIME STEP 

CHANGE 

CRANK-NICOLSON 45 CENTRAL 0.1 1.1 

1.2 

1.5 

0.01 1.2 
0.001 
0.0001 

UNIFORM 0.1 

SKEWED 

72 CENTRAL 0.1 1.0 

1.1 

1.2 

1.3 

1.4 



UNIFORM 

SKEWED 

132 CENTRAL 0.1 

UNIFORM 

SKEWED 

DVOGER 45 CENTRAL 0.1 

MTH=0 

0.1 



1.5 

1.2 

1.1 

1.2 

1.5 

1.2 



N. A. 



DVOGER 
MTH=1 
YMAX= 1 
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TABLE V (Continued) 

run method nodes disturbance convergence time step 

CHANGE 



24 








25 








26 


DVOGER 

MTH=1 

YMAX-1 


45 


CENTRAL 


27 






UNIFORM 


28 






SKEWED 


29 




72 


CENTRAL 


30 






UNIFORM 


31 






SKEWED 


32 




132 


CENTRAL 


33 






UNIFORM 


34 






SKEWED 


35 


IMPLICIT 


GEAR 45 


CENTRAL 


36 








37 








38 








39 






UNIFORM 


40 






SKEWED 


41 




72 


CENTRAL 


42 






UNIFORM 


43 






SKEWED 


44 




132 


CENTRAL 


45 






UNIFORM 


46 






SKEWED 


47 




45 


CENTRAL 



0.01 

0.001 

0.0001 



0.1 



0.1 

0.01 

0.001 

0.0001 

0.1 



0.1 



0.1 



0.5 
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TABLE V (Continued) 



RUN 


METHOD 


NODES 


DISTURBANCE 


CONVERGENCE 


TIME STEP 
CHANGE 


48 








1.0 




49 


DVOGER 

mth=i 14 

YMAX-10 ±4 


45 


CENTRAL 


i— i 

o 





TABLE VI 



DATA DECK ARRANGEMENT 



CRANK-NICOLSON 

Title 

NUMEL ,NUPBP ,NUMSNP ,NFULEL 
List of outer boundary points 
MTH , MAX DER,N COUNT 

ZNU , F IS FAC , HBAR , EPSVAL , ERRVAL , AFUEL 
TO ,H ,TF ,HMIN ,HMAX 

V - neutron velocity 
D - diffusion length 

SGA - absorption cross-section 

SGF - fission cross-section 

ALPHA - reactivity temperature coefficient 

PSIIV - initial disturbance flux 

System nodal points with axial and radial position 
Element nodal connectivity 
Nodal neighbor connectivity 

DVOGER 

Title 

NUMEL , NUPBP , NUMSNP , NFULEL 
List of outer boundary points 
MTH ,MAXDER,NCOUNT 

ZNU , F I SFAC , HBAR , EPSVAL , ERRVAL , AFUEL 
TO ,H,TF ,HMIN ,HMAX 

V 
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TABLE VI (Continued) 



D 

SGA 

SGF 

ALPHA 

PSIIV 

System nodal points with radial and axial position 
Element nodal connectivity 



IMPLICIT GEAR 



NCORE 

List of core elements 



with uniform initial disturbance 



Title 

NUMEL , NBP , NUMSNP , NFULEL , NELROW ,MXEVAL , NCOUNT , NCMPK 
NDE , NL 



List of outer boundary points 

VELOCT , DSUBF , DSUBC , SGMAAF , SGMAAC , SGMAF , FISFAC ,HBAR 
ZNU , ALPHA , RMSEPS , AFUEL , TSTART , TEND , PINITV , SCALE 
HMIN ,HMAX 

MTH ,MAXDER,NPROB ,NTYPE , JPLOT ,NPLOT , I RUN 



Nodal neighbor connectivity 

System nodal points with radial and axial position 
Element nodal connectivity 
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Z cm. 




Figure 1. Reactor Model with 45 Nodes 
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Figure 2 . Reactor Model with 7 2. Nodes 
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Z cm. 




Figure 3. Reactor Model with 132 Nodes 
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Figure 4. Nodal Neighbor Connectivity 
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0 cm., R = 60 cm. , 

cm. a Z = / 20 cm 
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Figure 5. Time Dependent Neutron Flux for All DOF Using Crank-Nicolson Method 
for a Uniform Disturbance Throughout the Core. 
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Figure 6. Time Dependent Neutron Flux for All DOF Using Crank-Nicolson 
Method for a Central Disturbance (R = 0 cm. , Z = 0 cm.) 
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Figure 7. Time Dependent Neutron Flux for All DOF Using Crank-Nicolson Method 
for a Skewed Disturbance (R s 60 cm., Z = 0 cm.) 
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Figure 8. Plot Comparing Degrees of Freedom to Computer Processing Time 
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Figure 9. Effect of Time Step Change in Crank-Nicolson on Achieving Steady State. 
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Figure 11. Time Dependent Neutron Flux for All Methods with 45 DOF 
for a Skewed Disturbance (R = 60 cm., Z = 0 cm.). 
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Figure 12. Time Dependent Neutron Flux for All Methods with 45 DOF 
for a Uniform Disturbance Throughout Core. 
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Figure 13. Time Dependent Neutron Flux for All Methods with 72 DOF 
for a Central Disturbance (R = 0 cm., Z = 0 cm.). 
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Figure 14. Time Dependent Neutron Flux for All Methods with 72 DOF 
for a Skewed Disturbance (R = 60 cm., Z = 0 cm.). 
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Figure 15. Time Dependent Neutron Flux for All Methods with 72 DOF 
for a Uniform Disturbance Throughout Core. 
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Figure 16. Time Dependent Neutron Flux for Crank-Nicolson and Implicit Gear 

Methods with 132 DOF for Central Disturbance (R = 0 cm. , Z = 0 cm. 
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Figure 17. Time Dependent Neutron Flux for Crank-Nicolson and Implicit Gear 

Methods with 132 DOF for Skewed Disturbance (R = 60 cm., Z = 0 cm. 
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Figure 18. Time Dependent Neutron Flux for Crank-Nicolson and Implicit 
Gear Methods with 132 DOF for a Uniform Disturbance 
Throughout the Core. 
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Figure 19. Comparison of Error Criterion for Implicit Gear Method. 
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Figure 20. Typical Approach to Steady State Solution Value 
with Various Error Criterion (Implicit Gear). 
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Figure 21. Time Dependent Neutron Flux for All DOF Using Implicit Gear 
Method for a Central Disturbance (R - 0 cm., Z = 0 cm.) 
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Figure 22. Time Dependent Neutron Flux for All DOF Using Implicit Gear 
Method for a Skewed Disturbance (R = 60 cm., Z = 0 cm.). 
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Figure 23. Time Dependent Neutron Flux for All DOF Using Implicit Gear 
Method for a Uniform Disturbance Throughout the Core. 
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Figure 24. Time Dependent Neutron Flux for 45 and 72 DOF Using the 
DVOGER (Gear) Method for a Central Disturbance 
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Figure 25. Time Dependent Neutron Flux for 45 and 72 DOF Using the 
DVOGER (Gear) Method for a Skewed Disturbance 
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Figure 26. Time Dependent Neutron Flux for 45 and 72 DOF Using the 

DVOGER (Gear) Method for a Uniform Disturbance Throughout 
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100 CONTINUE 
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WRITE ( 6 1 1 1 ) NUMEL, NBP ,NUM SNP , NFULEL , MELROW, MXE V AL , NCOUNT, NCMPK 
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WRITE (6 |40) 

40 FORMAT (/// 1< GEOMETRY OF SYSTEM NODAL POINTS'///) 
READ (5, 41) (SYSN|0D(I),R(I ),Z(I) , I=1,NUMSNP) 



FURMAT(5X,I5,2F15.5) 

WRITE (6,60) ( SYSNOD(I) ,R< I) ,Z( I ) , 1=1, NUMSNP) 

FORMAT (lX.'SVSTEM N3DAL POINT NO- ',I3,5X,'R COORD-= ',F10.6,5X 
IZ COORD. = ' , F 1 0 • 6 ) 
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STORE ELEMENT MATRIX, AMATRX , INTO SYSTEM MATRIX, BIGA 
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— 0ATA GENERATOR TO PROVIDE SYSTEM MODAL POINT 
RADIAL AMD AXIAL POSITION 

IMPLICIT REAL * 4 ( A-H , 0-Z) , I NT EGER * 4 ( I-N) 
INTEGER * 4 SfSNOD, ELNOD 

DIMENSION R( 1 5 3 ) , ZC 150) ,X(150) ,Y(150) ,SYSN0D(150) , 
1 ELNOD ( 220 » 3 ) 

READ IN GRID VALUES AND NJMBER OF NODES 

R EAD ( 5 1 100 ) NH,NV 
FORMAT (215) 

WRI TE ( 6 » 200) 

FORMAT ( ■ 1 • ) 

WRITE (6 t300) NH,NV 

FORMAT ( 2X t 1 NH =',I5,5X,»NV =*,I5) 

R EAD( 5 1 40 0) (X( I) ,I=1,NH) 

R EAD ( 5 ♦ 40 0) ( Y r ( J ) ,J=1,NV) 

FORMAT ( 15F5.1 ) 

WRITE (6 ,500) (X( I ),I=1,NH) 

WRITE (6 t501 ) (/( J) ,J = 1,NV) 

FORMAT ( //2Xt'R* , 5X , 14F8..3 ) 

FORMAT ( //2X, • Z*,5X,14F8.3) 

NVH=NV*NH 
NN=NV-1 
K = 0 

DO 10 1=1, NH 
DO 20 J = 1 , NN 
K =K+ 1 

SYSNOD { K ) =K 
R ( K) =X ( I) 

Z ( K) = Y( J ) 

CONTINJE 
CONTINUE 
K=NVH-NH 
DO 60 I =1 , NH 
K=K+1 

SYSNOD ( K ) =K 
R (K ) =X ( I) 

Z(K)=Y( NV) 

CONTINUE 
WRITE(6, 700 ) 

FORMA T( /////5 < , 1 NODE* ,1 IX , • R • , 1 IX , • Z • ) 

WRITE (6 ,600) (I , R( I ) , Z( I ) , 1=1, NVH ) 

WRITE (7,601) ( I , R ( I ),Z(I) ,1=1, NVH) 

FORMATt /5X, I 5,2F15.5) 

FORMA T(5X,I5,2F15.5) 



COMPUTE NODAL POINT CONNECTIVITY OF ELEMENTS 

LL= 1 
K = 1 
MM=0 

NHH=NH- 1 

N VV=N V- 2 

DO 30 1=1, NHH 

DO 40 J = 1 , NVV 

E LNOD (< , 1 ) =LL 

E LNOD ( K + l , 1 ) = wL 

ELN0D(K,2)=MMfNV 

ELN0D(K+1 ? 2)=LL+NV 

ELNOD (K,3)=ELN3D( K+l, 2) 

ELNOD (K+l, 3)=LL+1 

L L=LL+1 

K = K+2 

MM=MM+1 

CONTINJE 

LL=LL+1 

MM=MM+1 

CONTINUE 
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70 

750 

800 

950 



50 

900 

901 



00 70 1=1. NHH 
E LNOD (K , 1 ) =NN 
ELN0D(K+1 ,1)=NN 
ELN0D(K,2)=NN+NV-1 
E LNOD ( K+1,2)=LL+NV 
E LNOD( K,3)=ELN0D(K+1,2) 
ELN0D(K+1»3)=MM+NV 
K = K+2 

NN=NN+NV-1 

LL=LL+1 

MM=MM+1 

NEL=K-1 

CONTINUE 

WRITE(5 ,750) 

FORMAT! • 1 • ) 

WRITE ( 5 ,800) NEL 

FORMAT ( 5X , 1 NUMBER OF ELEMENTS =',I4) 

WR ITE (6,950) 

FORMAT ( /////10X, 'CONNECTIVITY MATRIX' ) 
DO 50 I = 1 » NEL 

WRITE(6, 900)1, (ELNOD(I ,J) ,J=1,3) 

WRITE (7, 901) I, (ELNDD(IfJ) ,J=1,3) 

CONTINUE 

FORMAT ( / 4 1 10 ) 

FORMAT (4110) 

IF(NV.EQ.12)ST3P 

GO T3 1 

END 



139 



ooo o ooo 



1 

100 

200 

400 

300 



10 

20 

30 

500 

600 



-NUCLEAR PROPERTY DATA GENERATOR- 



INTEGER * 4 NEL , ELEM, ITYPE 

DIMENSION V ( 220 ) , ALPHA ( 220 ) , SGF ( 220 ) , SGA ( 220) ,D(220) , 
1 I TYP E ( 220 ) » E LE M ( 220 ) 

R EAD ( 5, 100) NEL 

IDENTIFY CORE AND REFLECTOR ELEMENTS 
READ( 5,200) ( ITYPE ( I ) ,1 =1 t NEL) 

FORMAT (2110) 

FORMAT ( 1615) 

W RI TE ! 6 , 400 ) 

FORMAT ('l',2X,'ELEM',6X,'D',8X,» SGA' , 9X , 1 SGF • , 8 X , 

1 'ALPHA' ,8X, ' V ' ) 

FORMAT ( 2X, 14,5 ( IP E 12. 4) ) 

ASSIGN PROPERTY VALUES TO ELEMENTS 

DO 30 1=1, NEL 

ELEM( I) =1 

V( I ) = 4 • 8E+07 

ALPHAU ) =1 . 0E-05 

I F ( I TYP E ( I) • EQ. 0 ) GO TO 10 

D( I ) = 1 . 2 

SGA ( I ) = • 008 

SGF ( I ) = 0 

GO T3 20 

D ( I ) =0. 91 3 

SGA( I )=0.014 

SGF ( I ) =0 . 008 

WRITE! 6 ,300) ( I , 0 ( I ) , SGA ( I ) ,SGF ( I ) , AL PHA ( I ) , V ( I ) ) 
CONTINJE 

WRITE(7,500MV( I), 1 = 1, NEL) 

WRITE! 7 ,500) (0( 1 ) , 1=1 , NEL) 

WRITE ( 7,500) ( SGA! I ) ,1 = 1, NEL) 

WRITE(7 ,500) (SGF (I ), 1 = 1, NEL) 

WRITE ( 7,600) ( ALPHA! I ) ,1=1 ,NEL) 

FORM AT ( 8G1 0. 5 ) 

FORMAT ( 5F15 • 8 ) 

I F ( NEL • EQ . 64) 33 TO 1 
I F ( NEL . EQ . 1 12 ) GO TO 1 
I F ( NEL • EQ.220)ST0P 
END 



140 



ooor>o o oo onon 



■DATA GENERATOR TO PROVIDE NODAL \1E IGHBOR 
CONNECTIVITY 



INTEGER * 4 NV, NVH, LCON , NN, MM, LL , JJ , MNOD , I I ,KK, 
1 IJ,IK,IM,IN,NH,LI,LJ,LK 
DIMENSION MNOD( 132,7) 

READ IN INITIAL DATA 
I READ( 5, 200) N V, M VH , LCON , NH 
200 FORMAT (415) 

COMPUTE MATRIX COUNTERS 
M=N V+l 
N=NV-1 
N N=2*N 
MM=3*N 
I I =4*N 
KK=5*N 
I J=6*N 
I K=7*N 
I M = 8*N 
IN=9*N 
L K = 1 0 * N 
LI =NVH-NH 
L J=L I + 1 
LL=LI-NV+2 



CALCULATE COMPACTED NODAL NEIGHBOR CONNECTIVITY 
MATRIX 



10 



20 



30 



40 



50 



60 



DO 10 1=1. 
MNOD ( I , 1 ) = 
MNOD ( 1,2)* 
MNOD ( I , 3 ) = 
MNOD ( I , 4) = 
CONTINUE 
DO 20 1=2, 
MNOD( 1,5)* 
CONTINUE 
DO 30 I =NV 
MNOD( I , 6) = 
CONTINUE 
DO 40 I =M , 
MNOD( 1,7)* 
CONTINUE 
DO 50 1=1, 
MNOD ( I , 7 ) = 
CONTINUE 
DO 60 1=1, 
MNOD ( I , 6) = 
CONTINUE 
MNOD ( 1,5)= 
MNOD ( N, 2) = 
MNOD ( N , 4 ) = 
MNOD( NV , 5 ) 
MNOD ( NN , 2 ) 
MNOD ( NN , 4 ) 
MNOD( NN + 1 , 
MNOD ( NN +1 • 
MNOD ( MM ,2 ) 
MNOD ( MM , 4 ) 
MNODt MM + 1 , 
MNOD( MM+1 , 
MNODUI+l, 
MNOD ( NV H , 2 
MNOD ( II ,2) 
MNOD ( II ,4) 
MNOD ( 1 1 +1 , 
MNOD( II+l, 
I F( NVH. EQ • 
MNOD ( KK T 2 ) 



NVH 

I 

1 + 1 
N+I 
NV+I 

NVH 

1-1 

, NV N 
I-N 

NVH 

I-NV 

NV 

0 

N 

0 

0 

L I + 1 
L 1 + 2 
=0 

= L I +2 
= L I + 3 
5 ) =0 
7 ) =0 
= L I + 3 
=L I +4 
5 ) = 0 
7 ) = 3 
5 ) =0 
) =0 
=Lm 
= L I +5 
5 ) = 0 
7 ) =0 
45)30 
= L I + 5 



TO 15 
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15 

70 

80 

85 



95 



300 



90 

100 

101 



MNGD( K< ,4)=LIH 

MNOD ( KK + 1 » 5 ) = D 

MNOD( KK+1 * 7 ) =D 

MNOD ( IJ »2 1=11+5 

MNOD( I J 1 4 ) =L I +7 

MNOD (IJ+1»5)=D 

MN0D(IJ+1 ? 7)=D 

MNOD ( I K * 2 ) =L I +7 

MNOD ( IK » 4 ) =L I + 8 

MNOD( I K + l » 5 ) =0 

MNOD ( I K+l » 7 ) =0 

IF(NVH.EQ.72)3D TO 15 

MNOD ( IM » 2 ) =L I +3 

MNOD ( I M » 4 ) =L I *9 

MNOD( I M + l f 5) =0 

MNOD ( I M + l « 7 ) = 3 

mnod(lKt2)=lh-io 

MN0D(LK,4)=LI+11 

MNOD ( L K+l t 5 ) -0 

MNOD ( LK+1 1 7 ) = 3 

MNOD ( IN » 2 J=L I *9 

MNOD( IN,4)=LI+13 

MNOD( IN+i,5)=0 

MNOD( I N +1 » 7 ) = 3 

DO 70 I =LJ * NVH 

MNCD( I » 6) =0 

MNOD ( I » 7 ) -0 

CONTI NJ E 

DO 80 I =LL » L I 

MNOD ( I f 3 ) =0 

CONTINUE 

DO 85 I =LL ♦ L J 

MNOD ( I » 4 ) =0 

CONTINUE 

DO 95 I =L J » NV H 

MNOD ( I » 3 ) =N 

IF ( I • EQ *NVH) G3 TO 95 

M NOD ( I+1»4)=MNDD( I » 3 ) 

N=N+NV-1 
CONTINUE 
MNODdI f 2 ) =NV H 
WR 1 7E ( 6 » 3 00 ) 

FORMA T< 5X ,» HE X AGONAL CONNECTI VI TY • , / 
DO 90 1=1? NVH 

WRITE (6 *100) ( MNOD( I , J ) , J = 1 , L CON ) 
WRITE<7 ,101) ( MNOD ( I » J ) ? J = l, LCON) 
CONTINJE 

F0RMAT(/9(2X?I3) ) 

FORMAT ( 9 ( 2X r 13) ) 

I F (NV • EQ. 12 ) S TOP 

GO TO 1 

END 



1 NODE') 
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